clear 
est clear

********************************************************************************
*	Programs
********************************************************************************

* program to calculate outcomes for Kling Liebmann Bounds 
capture program drop get_KLBounds
program define get_KLBounds, eclass 
	
	*** Syntax 
	syntax varlist(numeric) , rounds(numlist) thresholds(numlist) treatment(varlist numeric) attrition(string) 
	
	*** calculate the mean and sd for T and C in each round 
	qui foreach v of local varlist {
		local r = substr("`v'",length("`v'"),1)
		tempvar  `v'_m 
		gen 	``v'_m'  = .
		tempvar  `v'_sd
		gen		``v'_sd' = .
		su `v' 						if `treatment' == 0   [aw=w`r']
		replace ``v'_m' = r(mean) 	if `treatment' == 0   
		replace ``v'_sd' = r(sd) 	if `treatment' == 0  
		su `v' 						if `treatment' == 1   [aw=w`r']
		replace ``v'_m' = r(mean) 	if `treatment' == 1   
		replace ``v'_sd' = r(sd) 	if `treatment' == 1  				
	}    
	
	*** Replace outcome vars 
	local t = 0 
	qui foreach thres of local thresholds {
	    local ++t
		foreach v of local varlist {
			gen 	`v'_kl_ub_`t' =  `v'
			replace `v'_kl_ub_`t' = ``v'_m' + `thres'*``v'_sd' if `treatment' == 1 & missing(`v')
			replace `v'_kl_ub_`t' = ``v'_m' - `thres'*``v'_sd' if `treatment' == 0 & missing(`v') 
			gen 	`v'_kl_lb_`t' =  `v'
			replace `v'_kl_lb_`t' = ``v'_m' - `thres'*``v'_sd' if `treatment' == 1 & missing(`v')  
			replace `v'_kl_lb_`t' = ``v'_m' + `thres'*``v'_sd' if `treatment' == 0 & missing(`v') 
		}
	}	
end

* program to calculate estimates for Lee Bounds 
capture program drop get_LeeEstimates
program define get_LeeEstimates, eclass
	*** Syntax 
	syntax varlist(numeric) , treatment(string) covars(varlist numeric) sample(string) tight(string)

	qui foreach v of local varlist {
	    *** get relevant round for sample indicator 
		local r = substr("`v'",length("`v'"),1)
		local sample_ `sample'`r'
		
	    *** Estimate 
		if ("`tight'" == "yes") 	 noisily leebounds `v' `treatment' if `sample_' == 1 [pw=w`r'], tight(`covars') 
		else if ("`tight'" == "no")  noisily leebounds `v' `treatment' if `sample_' == 1 [pw=w`r'], 
		
		*** Generate the matrices 
		*matrix list e(b)
		*matrix list e(V)
		matrix b1 = e(b)[1,1] , 0
		matrix rownames b1 = y1
		matrix colnames b1 = `treatment' _cons
		matrix b2 = e(b)[1,2] , 0
		matrix rownames b2 = y1
		matrix colnames b2 = `treatment' _cons
		*matrix list b1 
		*matrix list b2
		matrix V1 = e(V)[1,1] , 0 \ 0 , 0 
		matrix colnames V1 = `treatment' _cons
		matrix colnames V1 = `treatment' _cons
		matrix V2 = e(V)[2,2] , 0 \ 0 , 0 
		matrix colnames V2 = `treatment' _cons
		matrix colnames V2 = `treatment' _cons
		*matrix list V1
		*matrix list V2
		
		*** Create bogus estimations and replace matrices 
		* lower
		reg `v' `treatment', 
		erepost b=b1 V=V1 , rename
		eststo `v'_lee_lb
		
		reg `v' `treatment', 
		erepost b=b2 V=V2 , rename
		eststo `v'_lee_ub
	}
end 





********************************************************************************
*	I - Impacts of Cash + Therapy over just Therapy
********************************************************************************

est clear

use "${data}/Analysis_wide.dta", clear
distinct block club

*** vars in round 0 causing issue 
drop *0 *1

*** use weights, and do the fix 
drop w2 w3
* 10 CRs missing from sampling weights file... hence the assert(match master)
merge 1:1 cr_id using "${data}/PhaseII_status_weights.dta", assert(match master) nogen //keepusing(sample_weight)
tab sample_weight
*create a weight for each phase, set unweighted to 1 
gen w1 = 1  
gen w2 = 1 
label var w1 "Weights for RR, all equal to 1"
label var w2 "Weights for BL, all equal to 1"
* deal with w3 
gen w3 = 1 
replace w3 = sample_weight if !missing(sample_weight) 
* fix sample weight error 
replace w3 = 1 if w3 == 2
* remove unsampled folks 
replace w3 = . if phone_sample_PhaseII == 0 & field_sample_PhaseII == 0
label var w3 "Weights for EL, PhaseII selection probabilities"

*** drop control folks in longer follow ups 
drop if treat==0

*** create 'sample' tags 
gen sample2 = treat != 0
gen sample3 = treat != 0 & !(phone_sample_PhaseII == 0 & field_sample_PhaseII == 0)

*** attrition indicators 
forvalues i = 2/3 {
    gen att`i' =  missing(phq8_score`i')
	gen int`i' = !missing(phq8_score`i')
}

*** outcomes 
local varstoreg phq8_min? ghq12_min? phq8_score? ghq12_score? desired_fert? time_pref? paid_work? ed_enrolled? competencies? se_score? newlypregnant? newlymarried? riskysex?

*** Baseline estimates  
ds `varstoreg'
qui foreach v in `r(varlist)' {
	*get relevant round for sample indicator 
	local r = substr("`v'",length("`v'"),1)
	*center the blocks within round and non-missing sample
	capture drop blocks* 
	qui forvalues b=1/16 {
		gen blocks`b' = block == `b' 				if !missing(`v') & sample`r' == 1
		su blocks`b'				 				if !missing(`v') & sample`r' == 1 [aw=w`r']
		replace blocks`b' = blocks`b' - r(mean)		if !missing(`v') & sample`r' == 1
		su blocks`b'				 				if !missing(`v') & sample`r' == 1 [aw=w`r']
	}
	*regress 
	eststo `v' : reg `v' treat_cash b0.treat_cash##c.(blocks1-blocks15) if sample`r' == 1 [pw=w`r'], vce(cluster club)
}

 *** Get the KL bounded estimates 
get_KLBounds `varstoreg', rounds(2 3) thresholds(.1 .25) treatment(treat_cash) attrition("att")
ds `varstoreg'
qui foreach v of any `r(varlist)' {
	*get relevant round for sample indicator 
	local r = substr("`v'",length("`v'"),1)
	* centered blocks 
	cap drop blocks*
	qui tab block, gen(blocks)
	qui foreach b of varlist blocks* {
		su `b' 							if sample`r' == 1 [aw=w`r']
		replace `b' = `b' - r(mean) 	if sample`r' == 1
	}
	local blocks blocks1 blocks2 blocks3 blocks4 blocks5 blocks6 blocks7 blocks8 blocks9 blocks10 blocks11 blocks12 blocks13 blocks14 blocks15
	* regress
	eststo `v'_kl_lb_2 	 : reg `v'_kl_lb_2	treat_cash c.(`blocks')	b0.treat_cash#c.(`blocks') if sample`r' == 1 [pw=w`r'], vce(cluster club)
	eststo `v'_kl_lb_1 	 : reg `v'_kl_lb_1	treat_cash c.(`blocks')	b0.treat_cash#c.(`blocks') if sample`r' == 1 [pw=w`r'], vce(cluster club)
	eststo `v'_kl_ub_1 	 : reg `v'_kl_ub_1	treat_cash c.(`blocks')	b0.treat_cash#c.(`blocks') if sample`r' == 1 [pw=w`r'], vce(cluster club)
	eststo `v'_kl_ub_2 	 : reg `v'_kl_ub_2	treat_cash c.(`blocks')	b0.treat_cash#c.(`blocks') if sample`r' == 1 [pw=w`r'], vce(cluster club)
} 

*** get the Lee bounded estimates 
drop blocks*
qui tab block, gen(blocks)  
get_LeeEstimates `varstoreg', treatment("treat_cash") covars(blocks1-blocks16) sample("sample") tight("yes") 

*** Tabulate mental health outcomes 
local outcomes phq8_min phq8_score ghq12_min ghq12_score 
foreach r in 2 3 {
	local i = 0
	foreach v of local outcomes {
		
		* define labels 
		if 		("`v'"=="phq8_min") 	local coeflab = "PHQ-8$\leq$4"
		else if ("`v'"=="phq8_score") 	local coeflab = "PHQ-8 score"
		else if ("`v'"=="ghq12_min") 	local coeflab = "GHQ-12$<$3"
		else if ("`v'"=="ghq12_score") 	local coeflab = "GHQ-12 score"
		
		local ++i 
		* first time round replace
		if (`i'==1) {
			#delimit ;
			esttab 	`v'`r'_kl_lb_2 	 
					`v'`r'_kl_lb_1 	 
					`v'`r'_lee_lb 	 
					`v'`r'
					`v'`r'_lee_ub 	
					`v'`r'_kl_ub_1 	 
					`v'`r'_kl_ub_2 
				using "${tables}/bounds-mh-cash-`r'.tex", replace  booktabs
				label nonote nomtitle nonumber frag nolines nogaps noobs
				b(3) se(3) keep(treat_cash) coeflabel(treat_cash "`coeflab'")
				star(* 0.10 ** 0.05 *** 0.01)
				; 
			#delimit cr 
		}
		* then append
		else if (`i'>1) {
			#delimit ;
			esttab 	`v'`r'_kl_lb_2 	 
					`v'`r'_kl_lb_1 	 
					`v'`r'_lee_lb 	 
					`v'`r'
					`v'`r'_lee_ub 	
					`v'`r'_kl_ub_1 	 
					`v'`r'_kl_ub_2 
				using "${tables}/bounds-mh-cash-`r'.tex", append  booktabs
				label nonote nomtitle nonumber frag nolines nogaps noobs
				b(3) se(3) keep(treat_cash) coeflabel(treat_cash "`coeflab'") 
				star(* 0.10 ** 0.05 *** 0.01)
				; 
			#delimit cr 
		}
	}
}

*** Tabulate human capital outcomes 
* mid and endline
local outcomes ed_enrolled competencies se_score newlypregnant newlymarried riskysex
foreach r in 2 3 {
	local i = 0
	foreach v of local outcomes {
		
		if (!inlist("`v'`r'","se_score1","competencies2")) {
		
		* define labels 
		if 		("`v'"=="ed_enrolled") 		local coeflab = "Enrolled in school"
		else if ("`v'"=="competencies") 	local coeflab = "Competencies"
		else if ("`v'"=="se_score") 		local coeflab = "Self-efficacy score"
		else if ("`v'"=="newlypregnant") 	local coeflab = "Pregnancies"
		else if ("`v'"=="newlymarried") 	local coeflab = "Marriages"
		else if ("`v'"=="riskysex") 		local coeflab = "Risky sex"
		
		local ++i 
		* first time round replace
		if (`i'==1) {
			#delimit ;
			esttab 	`v'`r'_kl_lb_2 	 
					`v'`r'_kl_lb_1 	 
					`v'`r'_lee_lb 	 
					`v'`r'
					`v'`r'_lee_ub 	
					`v'`r'_kl_ub_1 	 
					`v'`r'_kl_ub_2 
				using "${tables}/bounds-hc-cash-`r'.tex", replace  booktabs
				label nonote nomtitle nonumber frag nolines nogaps noobs
				b(3) se(3) keep(treat_cash) coeflabel(treat_cash "`coeflab'")
				star(* 0.10 ** 0.05 *** 0.01)
				; 
			#delimit cr 
		}
		* then append
		else if (`i'>1) {
			#delimit ;
			esttab 	`v'`r'_kl_lb_2 	 
					`v'`r'_kl_lb_1 	 
					`v'`r'_lee_lb 	 
					`v'`r'
					`v'`r'_lee_ub 	
					`v'`r'_kl_ub_1 	 
					`v'`r'_kl_ub_2 
				using "${tables}/bounds-hc-cash-`r'.tex", append  booktabs
				label nonote nomtitle nonumber frag nolines nogaps noobs
				b(3) se(3) keep(treat_cash) coeflabel(treat_cash "`coeflab'")
				star(* 0.10 ** 0.05 *** 0.01)
				; 
			#delimit cr 
		}
		}
	}
}
*/







